;Calculate the centroids. We always do this in a 5x5 aperture...
function calc_centroids, im, lenslet_pitch, m_pix, flux=flux, dclamp=dclamp, sz=sz
if not keyword_set(dclamp) then dclamp=5 ;;A default value that seems very low but seems to work!
sep = lenslet_pitch/m_pix
psize = (size(im))[1]
if (lenslet_pitch gt 0.24 and lenslet_pitch lt 0.26) then begin
    n_lenslets = 12
    offset = floor(sep/2)
    posns = fltarr(2,n_lenslets)
    xy0 = psize/2 - offset
    posns[*,0] = [xy0,         xy0 -   sep]
    posns[*,1] = [xy0 + sep,   xy0 -   sep]
    posns[*,2] = [xy0 - sep,   xy0 ]
    posns[*,3] = [xy0      ,   xy0 ]
    posns[*,4] = [xy0 + sep,   xy0 ]
    posns[*,5] = [xy0 + 2*sep, xy0 ]
    posns[*,6] = [xy0 - sep,   xy0 +   sep]
    posns[*,7] = [xy0      ,   xy0 +   sep]
    posns[*,8] = [xy0 + sep,   xy0 +   sep]
    posns[*,9] = [xy0 + 2*sep, xy0 +   sep]
    posns[*,10]= [xy0      ,   xy0 + 2*sep]
    posns[*,11]= [xy0 + sep,   xy0 + 2*sep]    
    ;ic_jdm, im, /noc
    ;for i=0,11 do oplot, [posns[0,i], posns[0,i]], [posns[1,i], posns[1,i]], psym=2
endif else if (lenslet_pitch gt 0.15 and lenslet_pitch lt 0.18) then begin
    n_lenslets = 32
    offset = floor(sep/2)
    posns = fltarr(2,n_lenslets)
    xy0 = psize/2 - offset
    posns[*,0] = [xy0- sep,   xy0 -   2*sep]
    posns[*,1] = [xy0      ,   xy0 -   2*sep]
    posns[*,2] = [xy0 + sep,   xy0 -   2*sep]
    posns[*,3] = [xy0 +2*sep,  xy0 -   2*sep]
    for i=0,5 do posns[*,4+i] = [xy0 + (i-2)*sep,  xy0 -  sep]
    for i=0,5 do posns[*,10+i] = [xy0 + (i-2)*sep,  xy0       ]
    for i=0,5 do posns[*,16+i] = [xy0 + (i-2)*sep,  xy0 +  sep]
    for i=0,5 do posns[*,22+i] = [xy0 + (i-2)*sep,  xy0 +  2*sep]
    posns[*,28] = [xy0- sep,   xy0 + 3*sep]
    posns[*,29] = [xy0      ,   xy0 + 3*sep]
    posns[*,30] = [xy0 + sep,   xy0 + 3*sep]
    posns[*,31] = [xy0 +2*sep,  xy0 + 3*sep]
    ;ic_jdm, im, /noc
    ;for i=0,31 do oplot, [posns[0,i], posns[0,i]], [posns[1,i], posns[1,i]], psym=2
endif else stop
retvect = fltarr(2,n_lenslets)
flux = fltarr(n_lenslets)
if not keyword_set(sz) then sz=7
make_2d, findgen(sz)-sz/2, findgen(sz)-sz/2, xx, yy
for i=0,n_lenslets-1 do begin
    ;;Only centroid on the central 5x5 pixels. !!! Should be an option
    ;;more carefully thought out !!!
    subim = im[posns[0,i]-sz/2:posns[0,i]+sz/2, posns[1,i]-sz/2:posns[1,i]+sz/2]
    flux[i] = total(subim) > dclamp ;;Implement denominator clamping...
    xcent = total(subim*xx)/flux[i]
    ycent = total(subim*yy)/flux[i]
    retvect[*,i] = [xcent,ycent]
endfor
return, retvect
end
